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Abstract 

A general method is presented to explicitly compute autocovariance functions for non-Poisson 
dichotomous noise based on renewal theory. The method is specialized to a random telegraph 
signal of Mittag-Leffler type. Analytical predictions are compared to Monte Carlo simulations. 
Non-Poisson dichotomous noise is non-stationary and standard spectral methods fail to describe it 
properly as they assume stationarity. 
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I. INTRODUCTION 



A complete theory of stochastic processes was formulated by A.N. Kolmogorov in 1933 
[in. It included the characterization of stochastic processes in terms of finite-dimensional 

n 

distribution functions [2]. 

Since the early XXth century, it had been clear that the theory of stochastic processes 
has important applications in electronics and communications, for the basic understanding 
of devices such as vacuum tubes and, later, solid state diodes and transistors as well as for 
the behavior of cable and wireless communication systems. In particular, it turned out that 

noise can be described in terms by equations like the following one 

+00 

/(<)= Yl nt-tk) (i) 

k=— 00 

where F(t) is the current produced by an electron reaching the anode of a vacuum tube 
at time t, the k-th electron arrives at a random time and the series (PQ) is supposed to 
converge. S.O. Rice gave an account of the early developments in a Bell Labs monograph 
published in 1944 Q. 

Rice considers the example of the so-called random telegraph noise (RTN) or random 
telegraph signal (RTS), a stationary dichotomous noise where a current or a voltage ran- 
domly flips between two levels ±a as shown in Fig. 1 and flip events follow a Poisson process. 
An early treatment of random telegraph signals can be found in Kenrick's paper of 1929 



,4 



4j. According to Rice, Kenrick was one of the first authors to use correlogram methods 
to compute the power spectrum of a signal X(t). This was done by first obtaining the 
autocorrelation function 

R xx (^t; t) = E[(X(t) - E[X(t)])(X(t + At) - E[X(t + At)])}, (2) 

where E[X] denotes the expected value of the random variable X. Based on results by 
Wiener, Khintchin and Cramer, it is possible to show that, for a stationary signal when 
Rxx only depends on the lag At, the power spectrum is the Fourier transform of the 
autocorrelation function 5), Q]: 

/+00 
exp(iujv)R X x(v)dv. (3) 
-00 

For RTS with a = 1 and unitary average waiting time (a.k.a. duration), one gets 

Rxx =exp(-2|At|) (4) 



a results which will be derived below for positive At, leading to a Lorentzian power spectrum 

= i^-r (5) 

Even if very simple, RTS has interesting applications in different fields. A complete 
survey of the literature on RTS and its applications is beyond the scope of the present 
paper. However, it is interesting to remark that RTS has been used to explain 1/f noise 
in electronic devices. Indeed, A. Mc Whorter showed that an appropriate superposition of 
Lorentzian spectra for RTSs of different average duration gives 1/f noise [7^8, 9|. 



More recently, P. Allegrini et al. studied non-Poisson dichotomous noise 10] . Their paper 



is particularly relevant here, as a particular case of non-Poisson dichotomous noise is studied 
below in some detail, namely a generalized RTS where waiting times between consecutive 
level changes follow a Mittag-Leffler distribution. This particular process never reaches a 
stationary state in finite time and, for this reason, standard spectral methods fail to properly 
represent its properties. 

Our paper is organized as follows. In section II, devoted to theory, the stationary case 
is briefly reviewed. Then some tools are introduced taken from renewal theory. These 
methods lead to a general formula for the autocovariance function of a generalized RTS. 
This is given at the beginning of section III where the formula is then specialized to the case 
of Mittag-Leffler dichotomous noise. A comparison between analytical results and Monte 
Carlo simulations concludes this section. Section IV contains a summary of results as well 
as some directions for future work. 



II. THEORY 



A. The stationary random telegraph signal 

The random telegraph signal is a simple compound renewal process. A random variable 
X(t) can assume two opposite values, say X(t) = ±a, it starts with value, say, X(t) = +a 
for t = to = and it suddenly changes value assuming the opposite one at time instants 
ti, . . . , t n , . . . with the differences Tj = U — t^_i independently and exponentially distributed 
with common intensity fx. The probability density function of these differences is ip(r) = 
/xexp(— jir). If N(t) denotes the number of switches up to time t, the probability P(N(t) = 
n) of having n jumps up to time t is given by the Poisson distribution: P(n, t) = P(N(t) = 
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FIG. 1: Graphical representation of a RTS. 



n) = exp(—fj J t)(fit) n /n\. The random telegraph signal is a time-homogeneous stochastic 
process, therefore the autocovariance function: 



does not depend on t, but only on the lag At. 
B. Essentials of renewal theory 

A way to generalize exponential waiting times and the Poisson process is provided by 
renewal theory. A renewal process is a particular instance of one-dimensional point process. 
Renewal events occur at consecutive times to — 0, ti, . . . , t n , . . . and the differences Tj = 
U — tj_i are independent and identically distributed random variables. They are called 
durations or waiting times. The counting process N(t) related to a renewal process is the 



C xx (At;t) = E[X(t)X(t + At)} 



(6) 



positive integer giving the number of renewals up to time t. P(n, t) denotes the probabi 
of observing N(t) = n. Renewal theory is discussed in great detail in the book by Feller 



ity 



and in two monograp 
physicists also exists 



is, one by Cox and the other by Cox and Isham 12| . An account for 



The probability distribution P(n,t) can be related to the probability density, ip{r), of 
waiting times. It is sufficient to observe that t n , the instant of the n-th renewal event is also 
a random walk (in this case the sum of i.i.d. positive random variables): 

n 

*„ = $>, (7) 

i=i 

and the probability of having n jumps up to time t coincides with the probability of t n < t 
{n jumps up to t n and no jumps from t n to t): 

p( ni t) = i**r n ](t), (8) 

where * represents the convolution operator, ip* n {r) is the n-fold convolution of the density 
■0(t) and ^(t) = 1 — j^%l){u)du is the complementary cumulative distribution function 
(a.k.a. survival function). A continuous time random walk (CTRW) with exponentially- 
distributed waiting times and Poisson-distributed counting process is Markovian, moreover 
it has stationary and independent increments and, therefore, belongs to the class of Levy 



processes 
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151 ]. The exponential distribution is the only continuous memoryless 
distribution. In this case, the counting probability distribution P(n, At; t) of observing 
n jumps from a generic instant t to instant t + At coincides with P(n, At), that is the 
probability of observing n jumps from t = to At. In the general case, this is not true and 
P(n,At;t) explicitly depends on t [12]. In order to obtain P(n, At;t), one needs to know 
the distribution of the forward recurrence time a.k.a residual life-time y, the time interval 
between t and the next renewal event. If g(y,t) denotes the probability density function of 
the residual life-time, one has 

pAt 

P(n,At;t) = [g*^*r {n - 1) ](^t;t)= / g{y; t)P{n - 1, At - y) dy , (9) 

because the probability of n jumps occurring from instant t to instant t + At is given by 
the probability of t + tt n -i) < t + At and t + t( n _i) — t + J2i=i r «- The probability density 
function of the residual life-time g(y; t) can be written in terms of the renewal density h(t). 
The renewal density is the time derivative of the average number of counts up to time t, 



the so-called renewal function, H(t), that is the average value of the random variable N(t): 
H(t) := E[iV(£)], and h(t) := dH(t)/dt. The renewal density is a measure of the activity 
of the renewal process. For the Poisson process and exponentially distributed waiting times 
h(t) is a constant and coincides with the inverse of the average waiting time or the average 
number of events per time unit, in other words, h(t) = [L. It can be shown that [12j 

h(t)=ip(t)+ [ h(t - u)i){u) du , (10) 
Jo 



and that the density g(y; t) is given by [12] 



g(y,t) = rP(t + y)+ I h(t-u)4>(u + y)dy. (11) 



o 



C. A renewal process of Mittag-Leffler type 

F. Mainardi and R. Gorenflo, together with one of the authors of the present paper have 



161 ] . It is characterized by the following 



studied the renewal process of Mittag-Leffler type 
survival function 

*(r) = £,(- A (12) 
where Ep(z) is the one-parameter Mittag-Leffler function 



W = Ef^Tl) (13) 

where T(x) is Euler gamma function. The Mittag-Leffler function is a legitimate survival 
function for z = —r 13 , r > 0, < (3 < 1 and coincides with the exponential for [3=1. 
The Mittag-Leffler survival function interpolates between a stretched exponential for small 
values of r and a power law for large values of r. In particular, for < r << 1 one has: 

Ep(-T?) ~ 1 - ^ ~ exp(-^/r(/3 + 1)), (14) 



and for r — > oo 



H ^mm (15) 

7T TP 
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III. RESULTS 



A. Autocovariance for a general RTS 

Let us now consider a generalized random telegraph signal with durations that do not nec- 
essarily follow the exponential distribution. As the signal has zero mean, the autocovariance 
function Cxxif) = E[X(t)X(t + At)] coincides with the autocorrelation function. Through- 
out this paper, autocovariance functions will not be normalized. The generalized RTS is 
not a time-homogeneous (stationary) random process and the autocovariance function does 
depend on the initial time of evaluation. In particular, one gets 

oo 

C xx (At; t) = E[X(t)X(t + At)} = a 2 ^(-l) n P(n, At; t). (16) 

n=0 

Equation (|T6l) is a consequence of the fact that, in the time interval between t and t + At, 
with probability P(n, At;t), there can be an arbitrary, but finite, integer number n of re- 
newal events where the process of amplitude a changes sign. As P(n, At; t) is available 
from equations (Q, (jTOl) and ffTTT) . in principle, one can compute the time dependent au- 
tocovariance from the knowledge of the distribution of waiting times. However, the use of 
convolutions is painful and some analytical progress is possible by means of the method 
of Laplace transforms. Given a sufficiently well-behaved (generalized) function f(t) with 
non- negative support, the Laplace transform is given by 

POO 

f(s)= / eM-st)f(t)dt. (17) 



(A generalized function is a distribution in the sense of Sobolev and Schwartz |17l|). 

Here, the method of the double Laplace transform described by Cox turns out to be very 
useful 



121 ] . Let us denote by s the variable of the Laplace transform with respect to At or 
y and by u the variable of the Laplace transform with respect to t. Then g(s; u), the double 
Laplace transform of the residual life-time density g(y; t) in (fTTlh is given by 

g{s,u) = . (18) 

u — s 

Now, one gets from equation that 

P(n,s;u) = l(s;ums)r^(s) = W g ) ~ + h{u)) ^jn-ifr^ (19) 

u — s 



The double inversion of the Laplace transforms in equation (Fl9|) may be as a formidable task 
as the direct calculation of convolutions. However, equation ffl9l) can be used to investigate 
the asymptotic behaviour of P(n, At; t) for small and large values of t by means of Tauberian- 



like theorems If the average waiting time (r) = E[r] is finite, then for t — > oo, as a 
consequence of the renewal theorem, one has that h(t) — > l/(r) and the renewal process 
behaves as a Poisson process [12]. This is not the case for the renewal process of Mittag- 
Leffler type where h(t) = t p ~ l /Y{^) (and H(t) = tP/T(l + (3)). In other words, all the 
RTSs following a renewal process with E[r] < oo reach a stationary state after an initial 
non- stationary transient; in this stationary state, the autocovariance does no longer depend 
on t. The RTS of Mittag-Leffler type never reaches such a state. 

B. Analytical result for the RTS of Mittag-Leffler type 

Let us consider the case t = 0. Equation (fT6|) simplifies to 

oo 

C xx (At; 0\/3) = a 2 ^(-l)-P(n, At). (20) 

n=0 

Using the Laplace transform of the Mittag-Leffler survival function (fl2|) s^~ l /{l + s 13 ) and 
of the corresponding density function 1/(1 + s 10 ) Q, Q, one gets: 

C xx ( S] 0\(3) = a 2 ^^; (21) 

this Laplace transform can be inverted to yield 

C xx (At; 0|/3) = a 2 E (3 (-2(Atf). (22) 

For (3 = 1 this formula reduces to 

C xx (At; 0|/3 = 1) = a 2 exp(-2At). (23) 

Indeed, equation (1231 can be directly derived for a RTS with intensity /i = 1. In this case, 
P{n, At) = exp(—At)(At) n /n\, therefore from equation (|20|) one gets 

C X x(At; 0|/3 = 1) = a 2 exp(-At) = a 2 exp(-2At). (24) 

n=0 

Thus, equation ff22l) coincides with the autocovariance function for the usual RTS when 
= 1. 
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C. Monte Carlo simulations 



The essential ingredient for Monte Carlo simulation of a RTS of Mittag-Leffler type 
is the generation of Mi ttag -Leffler deviates. An efficient transformation method has been 
thoroughly discussed in [19] based on the theory of geometric stable distributions 



3, 



21 



Mittag-Leffler distributed random numbers can be obtained using the following formula 



23|: 



Gg£r» «*>)* <*> 

where u and w are real random numbers uniformly distributed in (0, 1). For /3 = 1 eq. (I25j) 
gives exponentially distributed waiting times with 

Once a method for generating Mittag-Leffler deviates is available, the Monte Carlo sim- 
ulation of an uncoupled CTRWs is straightforward. It is sufficient to generate a sequence of 
n(t) + 1 independent and identically distributed waiting times Tp t i until their sum is greater 
than t. Then the last waiting time can be discarded and n(t) alternate signs are generated. 
A Monte Carlo realisation starts, say, at X(0) = +a, then at time t\ = Tp t \ the process flips 
to the new value X(ti) = —a and stays there for t\ < t < t 2 with t% = r^i + rg )2 , and so on. 

In Fig. 2, Monte Carlo estimates of the autocovariance function Cxx{At;0\f3) are com- 
pared to formula (1221) for several values of j3 and a = 2. The estimates have been produced 
using an ensemble average estimator. For each of 10000 values of At, 1000 values of the 
product X(0)X(At) have been produced. The sampling frequency is assumed to be 1/1000, 
so that the total length of the signal is 10 arbitrary units (seconds). A direct eye inspection 
shows that the agreement between simulations and theory is excellent. 

Fig. 3 shows the effect of non-stationarity. The autocovariance functions Cxx{At;t\(3) 
are compared for {3 = 0.8, a = 2 and several values of t. For t = 0, the theoretical line is 
plotted. The dependence of the autocovariance function on the choice of t points to the fact 
that standard spectral methods routinely used in signal analysis would fail in this case. 



IV. DISCUSSION AND CONCLUSIONS 

The main result of this paper is given in equation ffl6l) for the autocovariance function of 
a generalized RTS, the double Laplace transform of the probability distribution P(n, At; t) 
is given in (Tl9l ; these results are further specialized to a RTS of Mittag-Leffler type in 
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FIG. 2: (Color online) Monte Carlo estimates of the autocovariance function Cxv(Ai; 0|/3) (red 
irregular lines) compared to the exact result in equation (|22p (blue smooth lines) for several values 
of (5 and a = 2. Time is in arbitrary units called seconds. The autocovariance functions are not 
normalized. 

equation (I22|) . The derivation of such equations is based on renewal theory. The analytical 
result in equation f|T6|) is then successfully compared to Monte Carlo simulations. 

An important feature of a RTS of Mittag-Leffler type is that the renewal density h(t) 
of the Mittag-Leffler renewal process vanishes for t — > oo for < (3 < 1. Therefore, the 
hypoteses of the renewal theorem are not satisfied and the stationary regime is never reached. 
Therefore, contrary to other renewal processes, such as the Weibull process, that, after a 
transient, become stationary, the renewal process of Mittag-Leffler type is non- stationary 
also if it is observed far from its beginning. In particular, standard spectral methods fail to 
describe the features of the RTS of Mittag-Leffler type. 

These considerations may be useful when studying the origin of 1/f noise, a problem not 
addressed here but which will be the subject of a future paper. Indeed, the renewal process 
of Mittag-Leffler type can be described as an infinite mixture of exponential distributions, 
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FIG. 3: (Color online) Non-normalized autocovariance functions C xx (At;t\(3) for (3 = 0.8, a = 2 
and several values of t. Time is measured in arbitrary units called seconds. 



for instance, one has for the survival function 

/>oo 

E P (~ tl3 ) = / M exp(-nt)dfM, 
Jo 

where 



#0) 



l 



sin(/37r) 



vr + 2 



Therefore, based on the ideas developed in 
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(26) 



(27) 



, one already expects to have a power 



spectrum following a power law. Using the fact that 

/+oo 
E p (-2\vf)exp(iuv)dv, 
-oo 

one gets 



(28) 



S xx (lu; 0\P) = 2Re[C xx (s = -ilu; 0\(3)], (29) 
and the power spectrum for t — can be evaluated from equation f[2~Tj) . It turns out that 

2 Alo 13 - 1 sin(/?7r/2) 



Sxx(u;0\P) = 2a z Re 



-iu 



2 + {-iuY 



4 + 4a;/ 3 cos(/57r/2) 



(30) 
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Once again, it must be stressed that, due to non-stationarity, the power spectrum depends 
on t and not only on the lag At. 
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